Subroutine thrmeq(rtemp,pres,phi,gtmp)

	USE ckvar
	
	INTEGER :: n,i
	REAL :: phi
	REAL :: rtemp,gtmp,pres,tmpmx,tmpmn
	REAL(dp) :: hreac,hprod	! Total enthalpy
	REAL :: totmol		! Total # of mole
	REAL :: errh
	tmpmx=4000.0
	tmpmn=600.0

	totmol=0.0
	do i=1,nspc
	totmol=totmol+nmspc(i)
	enddo
	mfi=nmspc/totmol	! initialize all speciesK	
	mfspc=mfi
	
!	write(*,*) '::::: Initial composition of reactants @phi=',phi,'::::::'
!	write(*,*) nmspc(1),spc(1),'+',nmspc(2),spc(2),'+',nmspc(3),spc(3)
!! end of reactant setup

!! Calculate total enthalpy of reactants
	CALL calthp(rtemp)
	hreac=0.0
	do i=1,nspc
	hreac=hreac+nmspc(i)*hk(i)
	enddo
!! end of calculation of enthalpy of reactant

	do 10 n=1,100
	gtmp=(tmpmx+tmpmn)/2.0
!! calculate Kp(T) for every reactions
	CALL calthp(gtmp)
	CALL calkp(gtmp)
!! Solve simultaneous equations using Jacobian Matrix !!
!	write(*,*) 'before, temp=',gtmp
!	mfspc=mfi
!	write(*,200) (mfspc(i),i=1,nspc)
	CALL calsimeq(mfspc,pres,phi,totmol)	
!	write(*,100) (spc(i),i=1,nspc)
	do i=1,nspc
	nmspc(i)=mfspc(i)*totmol
	enddo

!! Calculate total enthalpy of products
	hprod=0.0
	do i=1,nspc
	hprod=hprod+nmspc(i)*hk(i)
	enddo
	errh=abs(100.0*(hprod-hreac)/hprod)
!	write(*,*) 'errh=',errh
	if(errh.le.0.1) go to 20
	if(hprod.ge.hreac) then
	tmpmx=gtmp
	else
	tmpmn=gtmp
	endif
!	write(*,104) hreac,hprod,gtmp,tmpmx,tmpmn
   10	continue

   20	continue
	write(*,*) 'Adiabatic temperature=',gtmp
	write(*,*) 'errh=',errh,'Iteration#:',n
	return
  100	format(8(1x,a8,2x))
  104	format(5(1pe12.5,1x))
  200	format(8(1pe9.2,2x))	
end Subroutine
